Amorphous Systems in Athermal, Quasistatic Shear 
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We present results on a series of 2D atomistic computer simulations of amorphous systems sub- 
jected to simple shear in the athermal, quasistatic limit. The athermal quasistatic trajectories are 
shown to separate into smooth, reversible elastic branches which are intermittently broken by dis- 
crete catastrophic plastic events. The onset of a typical plastic event is studied with precision, and 
it is shown that the mode of the system which is responsible for the loss of stability has structure 
in real space which is consistent with a quadrupolar source acting on an elastic matrix. The plastic 
events themselves are shown to be composed of localized shear transformations which organize into 
lines of slip which span the length of the simulation cell, and a mechanism for the organization is 
discussed. Although within a single event there are strong spatial correlations in the deformation, 
we find little correlation from one event to the next, and these transient lines of slip are not to 
be confounded with the persistent regions of localized shear — so-called "shear bands" — found in 
related studies. The slip lines gives rise to particular scalings with system length of various measures 
of event size. Strikingly, data obtained using three differing interaction potentials can be brought 
into quantitative agreement after a simple rescaling, emphasizing the insensitivity of the emergent 
plastic behavior in these disordered systems to the precise details of the underlying interactions. The 
results should be relevant to understanding plastic deformation in systems such as metallic glasses 
well below their glass temperature, soft glassy systems (such as dense emulsions), or compressed 
granular materials. 



I. INTRODUCTION 

In crystalline materials it is generally accepted that the 
microstructural objects which govern deformation and 
flow are a class of topological defects known as dislo- 
cations. Most work in the field of crystalline plasticity 
focuses on describing deformation in terms of the under- 
lying dislocation dynamics. In the case of non-crystalline 
systems the situation is not so clear, even though a broad 
category of systems-including metallic glasses; clays and 
soils; pastes, foams, gels, and other so-called "soft" 
materials-seem to share a few hallmark traits. In the 
past several decades much work regarding the underlying 
microscopic processes of amorphous plastic flow, has left 
many unanswered questions and much controversy. Per- 
haps the most important of these questions is whether 
or not there exist some sort of microstructural defects in 
the materials which, roughly speaking, play the role of 
the dislocations in crystals P, S S S Q • 

Some of the earliest computer simulations of metal- 
lic glasses by Maeda and Takeuchi and co-workers Q 
and experiments on rafts of soap bubbles by Argon and 
Kuo ^ and Argon and Shi iQ] revealed that the atomic 
motions during plastic shear were confined to clusters 
of particles with a size of several particle diameters 
across. These observations inspired Argon to propose 
a theoretical scheme based on a mean-field treatment 
of transitions in local regions of space, "Shear Transfor- 
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mation Zones" (STZs), where each transformation con- 
tributes a quantum of plastic shear strain p^ . Draw- 
ing on Argon's works, Falk and Langer introduced 
mean-field equations of motion for the number density 
of STZs and showed that such a theory could account 
qualitatively for many aspects of the phenomenology of 
sheared metallic glasses such as strain hardening, the 
Bauschinger effect, and the emergence of a yield stress. 
The basic picture of an STZ may not be limited in 
utility to metallic glasses and could play an important 
part in the physics of other amorphous materials under 
shear such as foams, p ast es, gran ular materials and the 
like [ll|ll|ll[ll|li|ll|lllll|23 

These initial treatments all neglected spatial interac- 
tions of STZs. However, one might expect that the rear- 
rangement of an STZ should induce quadrupolar elastic 
displacements at long range, in an alogy with the trans- 
formation of an Eshelby inclusion, |2ll |22 | or the nucle- 
ation of a dislocation loop. Schematic, mesoscopic mod- 
els constructed to account for such elastic interactions, 
first proposed by Bulatov and Arg on, and later extended 
by others, [H ill HI El El El El have been shown 
to predict various sorts of localization of deformation. 
This mechanism might provide a physical explanation 
for an important technological problem: shear-banding, 
that is the localization of deformation to narrow bands, 
which is observed in many systems, including: metallic 
glasses 29, 30, 31], sheared rafts of bubbles 9], sheared 
foams confined between glass plates [s^, dry foams [ssf . 
and granular materials [SJ, IH Hg] . 

Despite these successes, theories of plasticity in amor- 
phous materials remain controversial because they rely 
on many assumptions which are difficult to check in 
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precise ways. In particular, these elusive STZs -unlike 
dislocations-cannot be identified a priori as any sort of 
topological defect ^6|. The athermal, quasi-static sim- 
ulations we present here provide us a starting point to 
perform such observations since they allow for the iso- 
lation and identification of elementary transitions be- 
tween mechanically stable states [s^, 113 • On these 
grounds, one might expect that an individual, elemen- 
tary, quadrupolar, Eshelby-like rearrangement might be 
associated with any particular single transition between 
mechanically stable states, but as we will see, the reality 
is more complex. Our work is a first attempt to simul- 
taneously identify elementary shear induced rearrange- 
ments in simulations of amorphous solids from both the 
perspective of the energy landscape and real-space. 

This paper is organized as follows. We first, in sec- 
tion m outline our athermal, quasistatic (AQS) algo- 
rithm, and provide a physical rationale for its use, dis- 
cussing the types of physical systems to which we expect 
our treatment to apply and its limitations. Section IIIII 
deals with the nature of the smooth elastic segments of 
AQS trajectories, and in particular how the trajectories 
break down at the onset of individual catastrophic events. 
The nucleation of one particular plastic event in a moder- 
ately sized system will be used as a case study. Section llVl 
deals with the nature of the individual plastic cascades 
themselves and their spatial structure. Again, a typical 
event will be used as a case study. Finally, in section Ivl 
we will show how the nature of the plastic events dis- 
cussed in section Hvl dictates particular scalings with the 
system size of the stress and energy relaxation during the 
plastic events and the strain interval between successive 
events. 



II. THE ATHERMAL, QUASISTATIC LIMIT 
A. Timescales 

Athermal, quasi-static (AQS) simulations have been 
used in several recent studies jH S HI El El El 
of plasticity in amorphous solids. The AQS algorithm 
simply consists of repeated alternating steps of: 1) min- 
imization of the potential energy of all the particles in 
the simulation cell [o^l and 2) application of a small, 
homogeneous strain to all particles and simulation cell 
boundaries. This simulation technique was introduced 
by Kobayashi, Maeda and Takeuchi 0, E3 as a way 
to bypass intrisic limitations of MD simulations to reach 
long timescales, and therefore low shear rates. Such lim- 
itations remain even with modern computers. 

The AQS algorithm relies on the idea that in the ab- 
sence of external drive, amorphous solids remain close 
to a mechanically stable state in a complex potential en- 
ergy landscape. For molecular or metallic glasses, this as- 
sumption is reasonable as soon as the bath temperature is 
low compared to the glass transition temperature Tg. A 
more precise bound can be obtained when these solids are 
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FIG. 1: A schematic representation of deformation-induced 
changes of a local minimum in the potential energy landscape. 
The shape of the landscape varies continuously as the strain 
is increased going from left to right with both the location 
and height of the minimum changing. 



submitted to some constant deformation rate 7, consid- 
ering that the thermal relaxation should be compared to 
7. The athermal limit correspond to the situation when 
7 » 1/Tr claxj where Tioiax characterizes the thermally 
activated escape of the system from a local minimum. 
In this limit, escape from local minima is primarily in- 
duced by strain and not by thermal activation [44L |45j . 
Because a low temperature limit is taken, this situation 
is likely to be relevant to different systems than metal- 
lic or molecular glasses. Foams, granular materials close 
to jamming, and several instances of soft glassy systems 
are intrinsically athermal and would likewise remain in a 
mechanically stable configuration in the absence of any 
external drive. 

When these amorphous solids are submitted to 
small amounts of deformation, they smoothly follow 
deformation-induced continuous changes of a local min- 
imum. This process is illustrated in figure ^ Strain in- 
duces a bias on the potential energy landscape, and the 
material configuration tracks the location of a single en- 
ergy minimum as it moves smoothly through configura- 
tion space. This kind of motion is completely reversible 
in that if we reverse the sense of the imposed strain, the 
system returns to its original configuration. As we will 
see, it is possible to solve analytically for the trajecto- 
ries of the system during this smooth motion, and these 
trajectories determine the elastic constants of the mate- 
rial. 46, 47, 48J Of course, reversibility holds only for 
small enough amounts of strain such that the minimum 
remains stable. For increasing strains, this smooth be- 
havior must eventually break down, as the energy mini- 
mum in which the system resides flattens out and collides 
with a saddle point HHHEIEI. 

From the preceeding picture, we understand that as a 
small strain rate is applied to such an athermal system, 
its response involves two types of behavior. Usually, the 
system smoothly and reversibly follows the continuous 
shear-induced changes of a single minimum as described 
above. Occasionaly, the occupied minimum vanishes, and 
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FIG. 2: Stress vs. strain curve for a 200x200 system of har- 
monic discs. The event at 7 = .1631 will be discussed further 
below in section HVI Note the smooth, roughly linear elastic 
segments interrupted by the discrete plastic events. 



the system has to relax toward an entirely new minimum 
in configuration space. This intermittent behavior shows 
up clearly when energy or stress is plotted against strain 
as in figure |21 On this figure, smooth segments corre- 
spond to reversible, elastic, changes of a particular energy 
minimum in configuration space; they are interrupted by 
discontinuous jumps, which corresponding to the shear 
induced annihilation of that minimum with a barrier. It 
is only during these jumps that energy is dissipated and 
across the jumps that irreversiblity may enter. 

The energy dissipation which occurs during these dis- 



crete events will take some finite time, 



ISSip. 5 



and 



order for the system to track the changes in the poten- 
tial energy landscape, it must be driven at a slow enough 
rate such that these discrete events have sufficient time 
to complete: 7 << l/rdissip.- Although the mechanisms 
of energy dissipation are system specific, it is reasonable 
to expect that material response in the quasi-static limit 
is largely determined by the existence of a potential en- 
ergy landscape and no t by the detailled mechanisms of 
energy dissipation. [5l|, 112, HI, |5J| ■ We will see, however, 
that in the AQS limit, the energy dissipation is strongly 
intermittent, and plastic jumps seen on figure |21 exhibit 
a broad range of amplitudes. As we will see, the typical 
amount of dissipation in an event, and accordingly the 
timescales which would be associated with that energy 
dissipation, show strong finite size effects. We should 
thus keep in mind that since Tdissip. might depend on the 
system size, it is probably only justified to speak of the 
quasi-static limit as a formal 7^0 limit, for a fixed 
finite system size. 

We see that the AQS limit entails three limits: zero 
temperature, zero strain limit, and the large size, thermo- 
dynamic limit. From the preceeding discussion it appears 
that the AQS limit holds when these limits are taken in 
the order: T ^ 0, then 7^0, then L ^ 00 



As long as the system remains in a convex region sur- 
rounding a minimum of the potential energy landscape, 
as it does along the continuous elastic branches, the pre- 
cise form of the energy minimization method should have 
no impact and the system will return to the local energy 
minimum after it is perturbed by the externally imposed 
deformation. On the other hand, when the system is 
driven past a limit of stability, as in the third frame of fig- 
ure ^ the precise method of energy minimization could, 
in principle, have an impact on the selection of a new 
minimum in which to reside as the system "rolls down- 
hill" away from the minimum which was just destroyed. 

The physical mechanisms for energy dissipation are 
modeled differently for different systems. In Durian's 
bubble model 55, 56], bubbles exert drag forces on each 
other proportional to their relative velocities; in Cundall 
and Strack's model for granular materials |57j . grains dis- 
sipate energy via a viscous dashpot connected in parallel 
with the springs which repel the particles; in simulations 
of molecular or metallic glasses jSg, one generally uses 
some sort of fictitious viscous thermostat to control the 
temperature in the system. Historically the AQS proce- 
dure has been implemented using some efficient energy 
minimization scheme, such as the non-linear conjugate 
gradient method, and we proceed along these lines. How- 
ever, one may hope that the details of the minimization 
technique, in particular, the nature of the viscosity and 
the existence of finite inertia, do not change the general 
picture that can be drawn from AQS simulations; a point 
of view we tentatively adopt here. 

This point of view finds some support in comparative 
studies of MD and AQS simulations, [H^l- Lacks has 
shown that the effective viscosity and diffusivity of the 
MD simulations extrapolates to the AQS results in the 
zero temperature limit for a low strain rate. In related 
work, Yamamoto and Onuki |3 have shown that 
the viscosity and diffusivity in similar simulations can be 
understood in terms of spatially heterogeneous dynam- 
ics which themselves are controlled by a critical point at 
T = 0,7 = 0. As we will see, these observation seem to 
be in agreement with ours and provide further indication 
that AQS simulation are a valid limit of MD simulations. 
To the best of our knowledge, however, no such explicit 
connection between MD and AQS has been shown for 
foam or granular models, but we consider it likely that 
analogous results would be obtained in these models in 
the limit of vanishing strain rate. We will therefore as- 
sume that the AQS procedure is equally applicable to the 
wet foams and frictionless granular systems, in addition 
to the metallic glasses for which it has traditionally been 
considered to apply. 



B. Numerical Details 

In this work we deal exclusively with 2D systems. Bi- 
disperse mixtures are used to inhibit crystallization. The 
mixtures used throughout have particles with radii: = 
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.5 and rs = .3 "9^ and a number ratio of: Nl = Ns^^^. 
The mixture is similar to that used by Falk and Langer 
and was found sufficient to inhibit crystalhzation. The 
systems are prepared via a zero temperature quench from 
an initiaUy random state as in references |^ |^. We 
suppose that the systems lose memory of their initial 
preparation before a strain of unity and use the strain of 
unity as a starting point for tabulating statistical prop- 
erties of the steady flow. No segregation is detectable 
within the 200% window over which the systems are 
strained, but it would be difficult to rule out effects which 
occur on very long strain scales. Lees-Edwards bound- 
ary conditions are used throughout ||5^]. The athermal 
quasistatic dynamics algorithm described above is imple- 
mented using strain steps of size 10~* which, as will be 
discussed below, is small enough such that all loading 
curves have well resolved elastic segments even for the 
largest samples simulated. 

Two different minimization algorithms from the conju- 
gate gradient family were employed [6^ : the non-linear 
PoUak-Riberi conjugate gradient minimizer and the trun- 
cated Newton linearized conjugate gradient minimizer. 
For the non-linear PoUak-Riberi conjugate gradient min- 
imizer, we used the routine as implemented in the GNU 
Science Library [6^. For the truncated Newton lin- 
earized conjugate gradient minimizer, we used an Armijo 
backtracking algorithm (621] with a sufficient decrease pa- 
rameter of .1 and the linear conjugate gradient routine 
as implemented in the Iterative Template Library [6^ 
adapted by us to perform truncation (62.] upon encounter- 
ing curvatures less than 10^^^. We found the latter pro- 
cedure to be more robust and efficient. The simulations 
of the single 50x50 system of Lennard- Jones particles dis- 
cussed in section Unl and the 200x200 system of harmoni- 
cally interacting particles discussed in section llVl utilized 
the former minimization algorithm, while the data runs 
which were used for the analysis in section^utilized the 
latter. Both algorithms gave statistically identical results 
when run on an ensemble of 50x50 harmonic systems. 

The potential energy functions were pairwise addi- 
tive central force laws. Three different force laws were 
employed: a standard 6-12 Lennard- Jones interaction, 
a harmonic, repulsive spring force, and a non-linear 
hertzian repulsive spring force (lOQj . The Lennard- Jones 
energy was truncated at a distance of 2 particle diame- 
ters, and linear and quadratic terms were added to ensure 
continuity up through second derivatives at the cutoff to 
avoid pathologies in the minimization routine. 

The pair interactions read: 

C^Harm(fii) = {I ~ Sij)'^ 9[1 ~ Si^) 

C/L,i(r.,) = (s,:^.i^-2sr6 + As.,+S4)0(2-,s,,) 

where is the unit step, and A and B are the coefficients 
which force continuity of the first and second derivatives 
at the cutoff in the Lennard- Jones potential. is the 



dimensionless separation between particles i and j, 
where Ri is the radius of particle i. 



III. ELASTIC BREAKDOWN AND PLASTIC 
NUCLEATION 

As described above, trajectories in AQS are composed 
of smooth, reversible elastic segments which are sepa- 
rated by discontinuous, irreversible plastic events. The 
smoothness of the elastic segments allows us to obtain 
analytical results regarding the singularity at the onset 
of a plastic event and to analyze the real-space structure 
of the critical mode of the system which is responsible 
for nucleating the subsequent plastic event. 



A. Analytical Framework 

We first briefly review the formalism developed in 
and recall the key results. The basic principle underlying 
this framework is that, owing to the smoothness of the 
energy during the elastic segments, exact analytical ex- 
pressions for the trajectory can be written by requiring 
that the system track the local energy minimum as its lo- 
cation in configuration space changes due to the imposed 
shear strain. 

The resulting equations governing the trajectory of the 
particles were found to be |lOlj] : 

""""^^^ -^»,73^^73 (1) 

The open circle over the r in equation is to indicate 
that the derivative is to be taken in the frame which is 
co-moving with the simple shear; that is: 

Hi) - Hi) (2) 

yii) = yii) - i^ii = 0). (3) 

Thus, equation describes the non-affine component 
of the motion, and the full motion of the particles in 
the laboratory frame is given by the imposed homoge- 
neous shear plus the correction term embodied in equa- 
tion . Hiajp is the so-called "hessian matrix" or "dy- 
namical matrix" - the second derivatives of the energy: 
Hiaj/j = Q^,^ ■ Sic is the derivative of the net force 
on particle i with respect to strain — or, equivalently, 
the derivative of the stress contribution of particle i with 
respect to a change in its position: = o^dr ■ 
cussed in detail in .48i], vanishes for configurations 
with local symmetry, and, as such, is a measure of the 
local configurational disorder about particle i. 
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The stress is defined as the total derivative of the en- 
ergy with respect to 7 while enforcing mechanical equi- 
librium of the particles via equation ^ . 

^^dU _ dU df„ ^dU _dU 
dj df-ia dj d"f d'f 

where the final equality holds because of mechanical equi- 
librium. So we see that the fact that the particles do not 
follow afHne trajectories does not give rise to any cor- 
rections to the stress. This is not true, however, for the 
shear modulus. 

To find the shear modulus, we simply differentiate yet 
one more time, again with a total derivative which should 
be understood to be taken while enforcing the corrections 
given in equation JQ). 

In equation lO, /ia is the term which arises from 

the Born expression for pure affine deformation |65l | , and 
/i„a = '^iaH~^^p'Ejf3 IS the term which comes from consid- 
ering the non- affine corrections. Note that the corrected 
modulus is always less than the naive Born expression. 
The microscopic analytical expressions for cr, /ia, Hiajp 
and Sict in the case of pairwise interacting systems can 
be found in reference 48], but we emphasize that the 
equations (Q), Q, and |SJ| are completely general and 
valid for any arbitrary n-body interaction potential (e.g. 
embedded atom methods |63 , or potentials for sili- 
con HI). 

So what is the structure in real space of the response 
during the continuous elastic segments? To illustrate, 
we now focus on a typical elastic segment in one par- 
ticular 50x50 Lennard-Jones sample. For any particular 
configuration, we can efficiently solve equation |^ via 
some iterative method. We use the conjugate gradient 
algorithm exa ctly as implemented in the Iterative Tem- 
plate Library [bJ with a relative tolerance of 10~*. This 
method, utilizing the analytical form for the elastic re- 
sponse, should be preferred over the alternate method of 
explicitly shearing the system by a small, finite amount 
then reminimizing the energy. The latter method es- 
sentially amounts to using a finite difference in lieu of 
a derivative whose analytical form we know. In refer- 
ence |43| , it was found using the alternate method of ex- 
plicit shear that when taking a small enough strain step 
to remain in the linear regime, the energy had to be com- 
puted to quadruple precision. No such special measures 
should be necessary in directly solving equation (^). 

Figure 13 shows the fields and Via at a value of the 
strain which is roughly a distance of 10~^ from the next 
catastrophic event. Note that the S field appears essen- 
tially random, as expected, based on its role as a measure 
of local configurational disorder. Via, on the other hand, 
should depend strongly on the low modes in the spec- 
trum of H, as can be seen directly from equation 




FIG. 3: The particular force in response to homogeneous 
shear, H (top), and the non-affine velocity (or "displacement") 
field, Via (bottom), at a strain configuration, 7 = 0.2945, or 
7, -7~ lO-''. 



It exhibits striking correlations in both compression and 
shear. 

The strongly spatially correlated behavior apparent in 
this elastic response field (and accordingly |/ino|) should 
be contrasted with the lack of any correlation beyond a 
length of a few particle diameters in any of the other 
mechanical quantities such as /Za, energy, pressure, shear 
stress, or vonMises stress. The importance of such non- 
affine corrections to elastic behavior has been realized 
in the context of experiments and simulations on emul- 
sions and foams reported by Liu and co-workers 69] and 
Langer and Liu [73 ■ More recently, Wittmer, Tanguy, 
and co-workers et. al. have conducted a comprehensive 
study of the contribution to elasticity of the non-affine re- 
arrangements in a Lennard-Jones system |46l |47| . How- 
ever, in this work, we will be more interested in the role of 
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FIG. 4: Stress (top) and shear modulus (bottom) for a small 
strain interval. Left: fixed strain steps of size 10~* using the 
standard energy minimization algorithm; Right: convergence 
to the yield point using a decreasing strain step size and the 
modified linesearch algorithm described in the appendix. 



We now proceed to isolate the singular behavior at the 
catastrophic points. In the lowest order non-trivial Tay- 
lor expansion for the energy at the transition point, we 
must include a higher order term for the critical direction, 
as the quadratic term vanishes. Generically, we expect a 
cubic term to remain [7l| : 

U ^ Asl + ^S.ia'Hiajl3Sji3+S'^{j: + EiaSia) + ^S'y^ (6) 

A,J^,T-liaji3,'^ia, and /Xq are constants which are evalu- 
ated at the transition. The physical interpretation of all 
but A are discussed at length in reference Sia de- 

notes the displacement of the reference coordinate away 
from the critical configuration: Sia = ridj) — ria{lc)- 
So denotes the projection onto the critical mode, so = 
Sia'^^^, where "^^^ is the unit vector which lies in the 
(non-translational) null space ofTCiajp. So we have: 

- = J. = 3^5^*0^ + HiajfiSjp + SjEia (7) 



the non-affine response on approach to the end-points of 
the elastic segments and its role in nucleating the plastic 
cascades rather than its role in renormalizing the elastic 
moduli away from the plastic events. 



B. Approaching Catastrophes 

Equation QJ, which describes the elastic segments, 
breaks down precisely when the local minimum vanishes; 
i.e. when the potential energy surface develops a direc- 
tion of zero curvature, along which the minimum collides 
with a first order saddle as in the cartoon of figure^ This 
scenario, with a single control parameter destabilizing a 
single degree of freedom, is the simplest possible type of 
bifurcation — known as a fold in bifurcation and catas- 
trophe theor y or a tangent bifurcation in dynamical sys- 
tems theory |7l| — , and we will see how this breakdown 
dictates the scalings of various quantities with strain at 
the onset of the plastic events. 

Since S depends only on local information about the 
near-neighbor particle configurations, we expect that it 
may be treated as roughly constant in the neighborhood 
of one of these catastrophic events, whence, according 
to equation the elastic response field will start to di- 
verge along the direction of the zero curvature. As the 
low curvature direction gets flatter and flatter, eventually 
the non-affine correction to the shear modulus dominates 
the Born term in equation ISJ and the net modulus be- 
comes negative. At this point the stress starts to decrease 
as a function of strain, and the system is unstable against 
any applied stress. These configurations are accessible to 
us because strain - and not stress - is controlled. Even- 
tually, the curvature will go to zero, and v and fina will 
diverge. 



USiaUSjp 

Requiring equation Q to be zero (which is equivalent 
to applying equation |^), and since lies in the null 
space of TLiajp, we have: 

Ti-iajpSip = (9) 

and 

iAsl = -S-fEo (10) 

We may solve equation Q for Sia up to the d uni- 
form translational modes and the critical mode. Equa- 
tion then provides the solution for the critical mode: 
So — \/— (^tSo/SA. The motion along this critical mode 
experiences a square root singularity, characteristic of 
the simple fold catastrophe, while the motion along the 
higher modes is smooth at the level of the expansion 

This singular behavior for sq induces singular behavior 
in the modulus and stress. Expanding equation (@J), we 
have: 



a ^ S - <57(S,„H-;-^S,;3) - + t^Jl (11) 

where S^q, is Sj^ with the critical component projected 
out: = SiQ, - ^'°^*°^Sj/3. Expanding equation 10), 
or, equivalently, taking the 7 derivative of equation l|ll|l. 
we get: 

M=^=Ma-Mna-(-^7r/^/|| (12) 

where /ia is the Born contribution to the modulus as 
in equation jSJl, fina = '^icJ^Jajp^jP is the non-singular 
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FIG. 5: All data shown here for the same trajectories plot- 
ted above in figure 2t>- a) relative participation of the low- 
est normal mode in the non-afhne elastic displacement field, 
a* = {"^laViaY / {viaVia) (dotted); lowBst eigenvalue of the 
dynamical matrix (solid); next several eigenvalues (dashed), 
b) In log-log scale (as a guide to the eye, the thick black line 
is y/^c — 7): (circles); lowest eigenvalue (squares); next 
several eigenvalues minus their terminal values (diamonds). 



part of the non-afRne correction, and the remainder is 
the isolated singular piece. 

To check these predictions, we perform a careful con- 
vergence to a particular catastrophic point in the system 
shown above in figure U] During the hnesearch portion 
of the minimization routine, we converged first to a min- 
imum of force and subsequently to a minimum of energy 
along the line and found this procedure to produce robust 
convergence to the transition point. 

Figure 0] shows the stress and modulus for the same 
system as in figure O upon approach to the singularity. 
Zooming in on the endpoint of the elastic segment, we 
see that the qualitative predictions of the theoretical ar- 
guments are borne out; namely that the stress reaches 



a maximum precisely as the modulus becomes negative. 
We further note that the Born contribution to the modu- 
lus is essentially constant on this region of 7 (not shown). 

In figure we plot the several smallest eigenvalues 
and the participation of the lowest mode, 0(7). In agree- 
ment with Malandro and Lacks, we find that a single 
eigenvalue vanishes. In the window of strain shown, the 
participation of the lowest mode goes linearly from about 
.85 at a distance of about 5^ ~ 10^"* to nearly unity at 
a distance of J7 ~ 10~^°. 

FigureEJ) shows the stress, the inverse of the modulus, 
and several higher eigenvalues as functions of strain on a 
log-log scale. The critical strain is determined to within 
10~^^ (not shown), and this terminal value is used to 
measure (57 = 7c — 7 for configurations down to 10~^° 
(shown). All quantities exhibit the same behavior 
at small 5^. This was predicted above for the modulus 
and critical curvature, however, the higher curvatures are 
constants at the order of the expansion ©. We can ra- 
tionalize the v^77 behavior for the higher modes by con- 
sidering that the total derivative of any function, f{sia), 
should be dominated by the singular behavior of sq: 



-1/2 



C. The Critical Mode 

The real space structure of a localized plastic event is a 
ke y inp ut into coarse-grained models of plasticity 22, 2^ 
milallilllllil and different supposed forms could lead 
to different emergent behavior in these models. A real 
space analysis of the incipient failure mode, although it 
does not correspond to a complete shear transformation, 
but rather to the onset of one, gives some insight into 
the form of an elementary plastic event. 

In figure El we show the field Via at Sj ~ 10"^", at 
which point it is almost entirely aligned with the crit- 
ical mode. The geometry is predominantly quadrupo- 
lar; a core region with outward particle velocities along 
the tensile axis and inward particle velocities along the 
compressive axis. We stress that the eigenmodes them- 
selves have changed little in the window of strain from 
(57 ~ 10"'* to (57 ^ lO^^*' as we have approached the 
incipient failure event, and the change in Via comes al- 
most entirely from the change in the relative weightings 
of the modes, with the quadrupole in the lower right of 
figure becoming dominant at the critical strain. 

In constructing figure |^ we have performed two trans- 
formations. First, we have applied an inverse affine trans- 
formation to the locations of the particles and their dis- 
placement vectors. Via'. 



Vy Vy 

VX ^ IVy 
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It is natural to ask whether the eigenmode is analogous 
to the displacement field generated in a linear, isotropic, 
elastic medium by some disturbance at the core. To map 
the discrete vector field onto a continuous field, we divide 
the system into annuli of width 2 starting at a radius of 2 
outside of the nominal core. In each annulus, we project 
onto circular harmonics by summing: 



vg{r; n) = ^ vgje 

jeAir) 



inOj 



(where j indexes the particles in a given annulus ) and 
normalizing each term by the number of particles in the 
annulus. 

It is found through this decomposition and is obvi- 
ous from simple visual inspection of the field, that the 
n — 2 (quadrupole) contribution dominates and has a 
phase angle which roughly gives extension along y = x 
and compression along y = —x (there is a slight clock- 
wise departure which can be seen in figure while the 
tangential component is roughly aligned at y = 0, x = 0. 

If the material behaves as a linear homogeneous 
isotropic elastic solid outside of some core region, we 
would expect, recalling the form of the quadrupolar space 
of solutions of the 2D Navier-Lame equation, that (72(| : 



Vr{r-2) = — + ^ '— 



veir;2) 



2A ^{1- k)B 



(13) 
(14) 



FIG. 6: The radial and tangential projections of the non-afiine 
displacement field, Via, at 7c — 7 ~ 10"^" transformed back to 
the rectilinear frame and centered on the vector with largest 
magnitude as described in the text. The center corresponds 
to the quadrupolar pattern which is starting to develop in the 
lower right of figure |21 at a strain of 7c — 7 ^ 10"'*. A core 
region of radius 2 is excluded from the image (and subsequent 
analysis). Scale is arbitrary: the field is normalized to unity. 
At this strain. Via is essentially indistinguishable from the 
critical eigenmode, '^ia- 



In this frame, the system maps onto itself under a purely 
vertical or purely horizontal shift by the box length, a 
property one which would be lost in a Lees-Edwards cell 
at finite strain. Further, we have moved to a co-ordinate 
system, (r, 6) , which is centered on the particle with the 
largest vector; for simplicity we have taken the center of 
the core to be on the particle with largest v. Also note 
that we have considered separately the radial and tan- 
gential projections and excluded a core region of radius 
2. 



where k is the ratio of Lame constants: n = ttttA — 7 and 
the 9 dependence of the radial and azimuthal fields should 
be understood to have a relative phase of 45 degrees. In 
figure[71 we plot the magnitude of the quadrupolar sector 
for Vr(r) and vg{r). 

The radial field seems to be consistent with a pure 
behavior, while the azimuthal field becomes too noisy to 
determine whether it follows any particular power law. 
For a general quadrupole, one would expect a crossover 
to r^^ behavior at small r, and it is evident that this 
crossover length is small if the r^^ term is even present 
at all. Furthermore, we may attempt to extract a value 
of K = {vr — vg)/{vr + ve). For radii less than 10, the 
azimuthal field is not too noisy, and the extracted k fluc- 
tuates between .4 and .5 (not shown), .4 being roughly 
consistent with the value of k averaged over the elas- 
tic segments. At larger radii, the extracted k leaves the 
physically allowed regime, becoming larger than .5, which 
is not surprising given the increased noise in vg at these 
radii. The agreement between this catastrophic mode 
and an elastic quadrupole is very reasonable. 
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FIG. 7: The magnitudes of the quadrupolar (n = 2) projec- 
tions of the radial, Ar2 (black circles) , and azimuthal Ae2 (red 
squares), components of of the critical mode. Solid lines are 
r^^ and r~^. 



D. Discussion 

There is a subtle distinction between the incipient fail- 
ure mode, which we have just measured, and a complete 
shear transformation. Ideally, one would like to measure 
changes in the particle configurations after the system 
has been driven past stability and then completely re- 
laxed. In our atomistic system, as in the coarse-grained 
models j^^, i22| , localized plastic events rarely occur 
in an isolated way, and most often the incipient failure 
event, such as the one measured here, triggers a subse- 
quent plastic cascade. As we will see below, although it is 
possible to make some measurements of the elementary 
shear transformations which occur during the cascades, 
they cannot be given as precise a meaning as the mea- 
surements of the incipient failure modes such as the one 
measured in this section. 

Before leaving the subject of plastic onset, we should 
clarify the relationship of the picture put forward here 
to the earlier numerical studies of Srolovitz and co- 
workers [t^. In those athermal simulations of metallic 
glasses under shear, the authors reported on stress con- 
centrations which acted as catalysts for plasticity. These 
stress concentrations were shown to vanish upon unload- 
ing, and an analogy was made with the stress fields 
around the tips of nascent shear cracks. Although we 
will, in fact, report on a related mechanism below, we 
have not found such stress concentrations at the on- 
set of plastic events. In fact, the formalism presented 
above (and discussed at length in |^) makes it clear 
that the instability is a collective property of the sys- 
tem and should not necessarily be discernible by look- 
ing at strictly local quantities. The local Born moduli, 
and the so-called "site-symmetry" parameters discussed 
by Egami et. al. 1741 . are closely related to the field 
(see reference |43 for the details of this relationship) 



which was shown to be essentially random and well be- 
haved near the transition. Although it is possible that 
there may be stronger correlations between Eia and the 
normal modes in other systems, such as those studied 
by Srolovitz and co-workers, no such correlations were 
found in any of the three different interaction potentials 
studied here in 2D, and there is no fundamental need for 
the stress itself to be particularly large locally in order 
to nucleate a plastic event. 



IV. PLASTIC EVENTS 

In the previous section, we focused on the behavior of 
the elastic segments of the quasistatic trajectories, and in 
particular, the behavior on approach to their endpoints 
at the onset of the plastic events. We now discuss the 
processes at work during these plastic events themselves 
as the system searches the potential energy landscape in a 
fully non-linear way in search of a new inherent structure 
after it is driven past a threshold of stability. A single 
typical event in a large system will be used as a case 
study. 



A. The Cascade Mechanism 

Already in the seminal work ^] , Argon emphasized the 
analogy between a local shear transformation and the 
nucleation of a dislocation loop. He cautioned that the 
analogy should not be taken too literally; in a crystal, the 
barrier for the nucleation of a pair of dislocations is large 
compared to the subsequent Peierls barriers, so, once nu- 
cleated, the pair is essentially free to glide apart. In a 
disordered system, on the other hand, even if one could 
topologically identify dislocations, the concept would be 
of limited utility, as there are no directions of symmetry 
along which the pair might glide apart. 

However, the elastic-like fields which are expected to 
result in the surroundings of a local shear transformation 
should alter the probabilities of observing subsequent 
shear transformations in neighboring regions. These 
elastic-like fields are expected to have quadrupolar sym- 
metry (i.e. they represent elementary shear) and resem- 
ble the elastic fields associated with the nucleation of a 
dislocation pair or the transformation of an Eshelby in- 
clusion 21]. 

Kobayashi, Maeda and Takeuchi were able to observe 
such elastic-like displacements of the particles surround- 
ing the core of a single shear transformation in their com- 
puter simulations |a,l3- Bulatov and Argon, with this 
picture of elastically mediated interactions in mind, con- 
structed a stochastic model of plasticity by embedding 
potential shear transformation sites unifmormily in a 2D 
lattice [23 13 US- Within their model, such cascades 
did emerge to play a role analogous to that of dislocation 
glide. Since then, several others 113, 113, [13, "2^] have con- 
structed coarse-grained models along similar lines with 
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FIG. 8: Schematic representation of a banding scenario. 
Sample is in tension along y = and compression along j: = 0. 
Stage 1) Virgin material. Stage 2) a) transforms and increases 
shear stress along red lines, increasing probability for b) to 
transform. Stage 3) Should the transformation at a) induce 
a transformation at b), the two stress fields are roughly addi- 
tive. 



varying assumptions about the precise microscopic de- 
tails of the elastic consequences of a local plastic event. 

The cascade mechanism is illustrated schematically in 
the cartoon in figure |S1 In frame 1, we have a vir- 
gin material with several potential shear transformation 
sites indicated by open circles. As a macroscopic strain 
is applied, eventually, the system will become mechani- 
cally unstable (as described in detail in the previous sec- 
tion), and a shear transformation will be nucleated at, 
say, region A. This initial event will have an associated 
displacement field with a cos{29) symmetry, giving in- 
creased (decreased) shear stresses along the lines y = x 
and y — —x (along the lines y — Q and a; = 0). The 
potential shear transformation sites which lie away from 
the initially transformed site along the lines of increased 
stress, having had their local stress levels increased, may 
themselves become mechanically unstable even without 
further increment of the macroscopic strain. The process 
favors lines of slip generated along the lines of maximal 
stress (which for our geometry would be along the ver- 
tical and horizontal axes); an emergent behavior much 
like the glide of a pair of dislocations, yet without the 
presence of any topologically identifiable objects. 



B. A Typical Cascade 

To fix these ideas, we recall and elaborate on the single 
event which was discussed in reference |75j . The system 
under consideration is a 200x200 system of harmonically 
interacting particles, with particle mixtures and prepa- 
rations as described in the previous section. A segment 
of the stress vs. strain curve is shown in figure |2| 

We examine the plastic event which occurred at 7 = 
.1631, and presume that it is representative of the kinds 
of events which occur at steady state. At 7 = .1631, the 
system has just been strained past the edge of an elastic 
segment - the local energy minimum has coalesced with 
a barrier, leaving only an inflection point in its wake - 
and the system is poised to undergo a non-linear, plastic 
rearrangement upon energy minimization. 

The total potential energy and sum of the squares of 
the forces during the energy minimization are shown in 
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FIG. 9: Energy and sum of the squares of the forces during the 
conjugate gradient descent for the single event at 7 = .1631. 



figure |5| As we use a conjugate gradient minimization 
routine, we have no rigorous notion of time. However, 
in steepest descent dynamics, where dxia/dt — Fia, the 
time derivative of the energy is precisely the sum of the 
squares of the forces. We have implemented a steepest 
descent dynamics and have checked on a single event that 
the resulting cascade was similar to the result of the con- 
jugate gradient algorithm. This allows us to interpret 
the number of minimization timesteps as some measure 
of time. However, the steepest descent algorithm is by 
far too slow to be used systematically in this study, and 
we thus have to rely on conjugate gradient methods. 

As the non-linear conjugate gradient algorithm pro- 
gresses, during the first 100 or so minimization steps, the 
system behaves much as if it was making a small cor- 
rection during a purely elastic relaxation, operating in 
an essentially linear regime. The energy relaxes toward 
a plateau, and the forces also decrease (not visible on 
the scale of the figure) accounting for the non-affine lin- 
ear elastic corrections (discussed in detail in the previous 
section) necessary to return the system, roughly speak- 
ing, to where the minimum "ought to" be. The forces 
are very small, and the energy plateau is very flat. It 
is almost as if the system is at a mechanical equilibrium 
state... a "quasi-basin" . Soon, however, at about the 
200-th minimization step, the forces get large and the 
energy drops rapidly; the system exits this quasi-basin. 

The force curve is intermittent with clusters of sharp 
peaks separated by quiescent periods, however it is diffi- 
cult to make precise distinctions; peaks may overlap and 
quiescent periods may be disturbed by small rumblings 
down by an order of magnitude from the peaks. We will 
focus on the cluster of force peaks which occur before 
minimization step 1150 and replot this segment of the 
descent on a blown up scale in figure [TUI Note the rea- 
sonably well defined initial force peak around step 230 
followed by several subsequent smaller peaks up to step 
800. After step 800, there is a period of relative quies- 
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FIG. 10: Close-up of the data shown in figure during the 
first 1150 minimization steps. 



cence which which lasts through step 1100. 

The energy descent is broken into intervals of 50 steps 
each, starting at 200 and ending at 1150. The slip 
and energy dissipation in realspace which occurs during 
these periods is shown in the sequence of images in fig- 
ures ^2 ^1 E^nd ^1 For any particular particle, we de- 
fine the slip as the difference between the displacement of 
that particle and the average displacement of its neigh- 
bors (particles with which it is in contact); in this sense 
it is analogous to a discrete Laplace operation on the dis- 
placement field. All real space structures are transformed 
back to the rectilinear frame by an inverse simple shear of 
magnitude equal to the current strain, .1631, such that in 
the plots shown, the point {x, y) may be identified with 
the points {x + aL, y + bL) where a and b are arbitrary 
integers. 

The incremental slip occurring in each window of 50 
minimization steps is shown in figure ITTI and the cumu- 
lative slip (with incremental slip superimposed in red) 
is shown in figure ^1 Any particle is shaded if it has 
slipped by more than lO"'^. The local incremental energy 
dissipation is plotted, in figure [T51 on a log scale which 
ranges from 10~^ to 10^, with grey indicating an energy 
change of less than 10~^ in magnitude, white indicating 
an energy increase, and black, an energy decrease. 

The pattern which emerges is that each of the peaks in 
the squared force in figure [TUl corresponds to a well local- 
ized cluster of particles which undergoes large slipping 
relative to its neighbors. As can be seen in figures ITTI 
and 1 121 the new slippage does not occur precisely on top 
of the slippage which occurred during previous peaks, 
but instead, new slippage tends to occur at the extremi- 
ties of the region which has already undergone apprecia- 
ble slip. Furthermore, the local energy dissipation which 
occurs concomitantly with each force peak and cluster 
of slipping particles can be seen to take the form of a 
quadrupole which is centered over the cluster of slipping 
particles. Note that the quadrupoles predominantly have 
energy increases (white) along the tensile axis, y = x, 
and energy decrease (black) along the compressive axis, 
y = —X. This spatial organization can be understood in 
terms of the banding argument outlined above. 

Figure ^] shows the displacement and slip which oc- 
curs after the energy relaxation is fully complete. The 
displacement field is plotted such that an arrow of length 
.5 is drawn at each particle in the direction of its displace- 



ment. The shade of the arrows represents the amplitude 
of the displacement on a linear scale. This representa- 
tion allows for better appreciation of the orientation of 
the displacement field even when its magnitude becomes 
small. Recall from the definition above that the slip is es- 
sentially a type of discrete derivative of the displacement. 
We prefer to deal directly with displacements rather than 
energies or stresses, as the latter are essentially spatial 
derivatives of the former and thus much noisier. 

The coherent, elastic-like behavior of the displacement 
field is striking. It is consistent with the banding mech- 
anism described above, and is roughly equivalent to the 
displacement field resulting from the gliding apart of a 
pair of dislocations or a shear crack running vertically 
through the system. The cascade appears to arrest be- 
fore it has spanned the system. There is a large, weak 
vortex (clockwise) at about {x — 0,y = L/2) between 
the slip line and its periodic images and a large, weak 
hyperbolic flow at about {x — 0,y = 0). They are to 
be expected from a periodic array of incomplete vertical 
slip lines. If the slip line had been complete, the result- 
ing pattern would have likely been the text book example 
of a pure shear with displacements in the y direction, a 
gradient along the x direction, translationally invariant 
along the y direction, and with a discontinuity along the 
slip line. 

Along the vertical slip itself, one sees small vortex-like 
(counter-clockwise) displacement fields in which the mag- 
nitude of the displacements is markedly smaller than the 
displacements of particles which slip. These smaller vor- 
tices represent regions of the material along the slip line 
which have failed to slip and can be thought of as result- 
ing from a uniform line of elastic quadrupoles (a slip line) 
with a gap in the line. Not surprisingly, these vortex-like 
regions of "unbroken" material which appear in the dis- 
placement field in figure 1141 correspond to regions which 
are absent from the plot of the slip field in figure [T51 



C. Discussion 

In this section, we have examined one particular, typ- 
ical (the magnitude of stress relaxation was around the 
average, and the event was taken at random from the set 
of all plastic events) plastic event in a 200x200 system 
of harmonically interacting particles and shown that the 
behavior was consistent with the picture of a cascade of 
elastically interacting localized plastic yielding events or- 
ganized into a line of slip along the vertical bravais axis of 
the simulation cell. We have observed other such events 
aligned along either the vertical or horizontal axis of the 
simulation cell. Such organization of local slippage events 
into lines of slip was observed longago in experiments on 
bubble rafts by Argon and Kuo ^ and, more recently, 
in confined films by Abd el Kader and Earnshaw |33| . 
Although experimental observations of these features are 
limited to studies of soap bubbles, observation in numer- 
ical simulations are quite rich. 
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FIG. 11: Incremental slip (as defined in the text) at 50 step intervals during the minimization routine. Sequence starts at step 
200 and ends at step 1150. A particle is shaded if it slips by more than 10^''. 



In their vertex model for dry foams, Okuzono and 
Kawasaki jj^l observed that in the limit of small strain 
rate, their system underwent large events, and the au- 
thors proposed an analogy with avalanches in sandpile 



models. The displacement fields associated with these 
events were similar to the ones we have shown here, with 
the system exhibiting two "elastic-like" regions slipping 
with respect to each other along a line at 45 degrees to 
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FIG. 12: Cumulative slip (incremental slip superimposed in red) at the same intervals as in figure ITTI 



the principle axes of applied strain (c./. figure 5b in ref- 
erence [73 ) . We do not know of any attempts to observe 
such displacement fields in models of wet foams such as 
Durian's bubble model jH^ |^ . 

Evidence of such transient, slip bands is also observed 



in the simulations of compressed granular materials by 
Aharonov and Sparks |3| and Kuhn [s^ [s^. These 
simulations took into account the Coulombic friction be- 
tween particles at contact and were performed at finite, 
but small strain rates, yet the emergent behavior seems. 
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FIG. 13: Local energy dissipation during the same time intervals as in figures inland [El Data is on a log scale which ranges 
from 10~® to 10^, with grey indicating an energy change of less than 10~^ in magnitude, white indicating an energy increase, 
and black, an energy decrease. 




FIG. 14: Particle displacements which occur during the entire plastic event. Individual arrows have a uniform length of .5 and 
a shading which is linear in the amplitude of the displacement. The reader is encouraged to utilized the zooming features of 
the pdf document format to explore the fine scale structure of the displacement field. 
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FIG. 15: Local slip (as defined in the text) which occurs during the entire plastic event. Arrow lengths are equal to the 
magnitude of the particular slip scaled by a factor of 10. Slips of amplitude less than 10"'^ are not shown for clarity. 



at least qualitatively, insensitive to these details. Note, 
however, that one might expect qualitatively different be- 
havior as the truly rigid, hard sphere limit is approached. 



Simulations were performed by Leonforte and co- 
workers |77| which were similar to the ones we have de- 
scribed here but with rigid walls at y = and y = L 
rather than fully periodic boundary conditions. Not sur- 
prisingly, the vertical bands are suppressed, and transient 
horizontal bands emerge. 



V. FINITE SIZE SCALINGS 

The cascades discussed in the previous section which 
comprise the plastic events will now be shown to give rise 
to simple scalings of various measures of the event size 
with the length of the simulation cell. In the following, we 
classify each strain step taken in the simulation as either 
elastic or plastic: if the stress was positive, a strain step 
was considered to be plastic if it resulted in an energy de- 
crease; if the stress was negative, a step was considered to 
be plastic if it resulted in an energy increase 

[loaj; oth- 



erwise, it was considered to be elastic. We suppose that 
this operational definition corresponds closely to the rig- 
orous definition of a plastic event in terms of the onset 
of a catastrophe via the vanishing of an eigenvalue of the 
hessian matrix, but as we emphasized above, care must 
be taken to ensure that we take small enough strain steps 
such that we remain in the quasistatic regime. Note how- 
ever that, in principle, there may always be small plastic 
discontinuities whose energy drop would be masked by 
the elastic energy increase for a given strain step size. 
All the data used for the statistical analysis were gath- 
ered from an ensemble of 8 systems at strains between 1 
and 2 for each box size and interaction potential. 

Before we begin our discussion of the statistics of the 
plastic events, we first report, in figure [TBI on the distri- 
butions of the instantaneous shear modulus, fi = for 
the various system sizes and interaction potentials. We 
normalize the values by the average flow stress, (a) and 
note that the distributions are essentially unchanged if 
we include /exclude the plastic events from the analysis. 
The average flow stresses for all interactions and system 
sizes are reported in table For each potential, there is 
a slight increase with size of the average, (a), which is 
noted in table U — this effect seems to saturate for the 
larger systems. 

The distributions of the moduli are all essentially 
Gaussian in the neighborhood of the most likely value, 
but with tails on the soft side, indicative of the non- 
linear yielding upon approach to the catastrophes dis- 
cussed above. The peak occurs at a dimensionless mod- 
ulus of about 28 for the Harmonic and Lennard-Jones 
systems and at a slightly higher value for the Hertzian 
system. This implies that, generically, these dense, disor- 
dered systems have a flow stress which is about .035 times 
the characteristic modulus, regardless of the nature of the 
interaction potential. This value of the "flow strain" , , 
is roughly consistent with simulations of various models 
of foams in 2D [H HI [71 S Iz3 (in the dense hmit, 
away from the loss of rigidity) and the 3D simulations of 
yielding of a Lennard-Jones glass 1^ , although 
it is somewhat higher than Johnson's recently reported 
universal flow strain in 3D glasses js^ . It is not clear pre- 
cisely how this flow strain, measured in steady flow, re- 
lates to the yield strain measured by Mason and cowork- 
ers in_c«cillatory rheological experiments on microemul- 
sions 84] , but the order of magnitude of a few percent is 
in rough agreement. 

We stress that it is our ability to resolve the elastic 
behavior and measure a well deflned modulus which gives 
us confidence that we have chosen a strain step which is 
small enough to properly resolve the quasistatic behavior. 
For larger strain steps, the well defined peak disappears, 
the stress essentially makes a random increase or decrease 
at each step, and the quasistatic behavior — i.e. the 
separation of plastic from elastic events — is lost. 

For the purposes of this work, the more interesting 
quantities are the various measures of cascade size. Ob- 
vious choices are the distributions of the stress drops, the 
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TABLE L Average stress during steady flow for the three 
different system sizes and interaction potentials. 
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FIG. 16; Distribution of instantaneous modulus (defined by 
the finite difference, Aa/Ay) normalized by the average stress 
during steady flow. Each group of three curves corresponds 
to a particular system length (12.5,25,50) arbitrarily shifted 
vertically for clarity with increasingly longer lengths shifted 
upward. Legend is as follows: circles with solid fine (black): 
Harmonic; squares with dotted line (red): Hertzian; diamonds 
with dashed line (green): Lennard-Jones. 

energy drops, and the lengths of the elastic segments. 
As the patterns we observed in the 200x200 system of 
harmonic discs were clearly ID features, one might ex- 
pect that the distributions of the various quantities which 
characterize the event size would be invariant when prop- 
erly rescaled by the length of the box. 

First consider the distribution of stress drops. Since 
the events we observe are predominantly organized into 
lines of slip which extend across the length of the cell, 
we expect that such an event should release an amount 
of stress equal to Act ~ ii6e ^{a/L) where a is some 
measure of the amplitude of the slip at the site of the 
cascade and L is the length of the box. Since the system 
is in steady state, the stress built up in the strain interval, 
/iA7, between these events must be equal, on average, to 
the stress released in an event, Act, so we must have that, 
A7 ^ a/L. The stress drops can be related to the energy 
drops if we make the simple assumption that, on average, 
the energy release is elastic; that is: 

AU - L^— aL{a) 

A* 

Now that we have considered the energetics of the plas- 
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tic events, we show how the same size effects dictate 
their relative frequency. During the elastic segments, the 
stress increase is ActoI ^ M^7- I^i steady state, this elas- 
tic stress increase must be balanced, on average, by the 
plastic dissipation, so we simply have: 

Act a^A 

A7 r- 

So we see that the relative frequency of the plastic events 
should not even depend on the energy scale of the under- 
lying interaction potential. 

We now proceed to plot the various distributions: 
P(;^),P(Acr(^)), and P{LAj), for all three system 
sizes and all three interaction potentials. All 27 curves 
collapse onto a single master curve which is very well 
fit by an exponential (with better collapse for the larger 
system sizes). For clarity, in figure [T7I we first show the 
distribution for each of the three quantities — energy 
drop, stress drop, strain interval — in its own plot for 
all 9 systems, then, in figure ITSl show the 9 curves cor- 
responding to only the largest system size, but for all 
three quantities and for all three interaction potentials. 
This characteristic value for a which one extracts from 
the master curve is a few tenths of a particle diameter, 
in good agreement with the discontinuity in the displace- 
ment fields shown above. The precision of the collapse in 
figure 1181 is quite striking and indicates that the nature 
of the slip which occurs along the cascade line has little 
to do with the precise nature of the interaction potential. 

The event rate scalings, in particular, have profound 
consequences concerning both fundamental and technical 
issues. For any arbitrarily small strain step size, there 
will always be systems large enough, such that the step 
size is no longer small enough to resolve the elastic behav- 
ior, many independent plastic events will be simultane- 
ously nucleated at every strain step, and the quasistatic 
behavior will break down. Thus, it is technically impor- 
tant, for any quasistatic simulation, to verify that one is 
properly resolving the elastic behavior, independently for 
each system sizel For the present study, this is evidenced 
by the peak-like behavior in figure^l We note that, had 
the elementary relaxation events been uncorrelated, the 
plastic event rate would have scaled extensively, as L^, 
so the correlations reduce the frequency of plastic events 
relative to what would have been observed for indepen- 
dent uncorrelated events. 

VI. CONCLUSION 
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FIG. 17: Distribution of (top to bottom): -^§y, and 
I/A7. Where AU is the energy drop for a plastic event. Act 
is the stress drop for a plastic event, and A7 is the length 
of an elastic segment. Legend is as follows: solid (black): 
L = 12.5; dotted (red): L = 25; dashed (green): L = 50; 
circles: Harmonic; squares: Hertzian; diamonds: Lennard- 
Jones. 



In conclusions, we have shown that elementary plastic 
events take the from of cascades of shear transformations 
in the athermal quasistatic limit for a general class of 
densely packed simple amorphous materials. The incre- 
mental stress and displacement fields associated with the 
cascades are quite reminiscent of micro-structural shear 
cracks which immediately heal themselves after crack- 
ing. These local shear transformations might be roughly 



thought of as the elementary objects responsible for plas- 
ticity in amorphous material, analogous to the disloca- 
tions which are thought to be responsible for plasticity 
in crystals. 

The locus of rearrangement shows some spatial corre- 
lation from one event to the next, but these correlations 
decay quickly after just a few plastic events. We ob- 
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FIG. 18: Probability distributions of and LA7 

for all three interaction potentials. Only the 9 curves corre- 
sponding to the largest systems are shown for clarity. 



serve no evidence for the kinds of pronounced persistent 
shear localization which is seen in many experiments and 
simulations of sheared amorphous materials where hard 
walls are employed to drive the system, and we find it 
likely that the persistent localization observed elsewhere 
is due largely to elfects of the boundary. Future inves- 
tigations into the role of the boundaries will be crucial 
in both atomistic studies [Sfil Is^ Is^L Is^ and mesoscopic 
models 22, 87]. 

The detailed picture of the onset of individual plas- 
tic events which we developed showed that directions 
of low curvature on the potential energy surface, when 
driven to zero by the strain, are responsible for nucleat- 
ing the cascades of shear transformations. These low- 
lying modes which were shown to have an essentially 
"shear-transformation-like" quadrupolar character, were 
found to play a dominant role in the non-affine elastic dis- 
placement fields, even reasonably far away from the on- 
set of a cascade, and might be observed experimentally 
in this way. Strikingly, at least in the systems studied 
here, these modes are not observable via looking at local 
stresses (e.g. as in reference JSj); '^^ local Born values of 
the elastic constants, and can only be observed through 
the non-affine elastic response which is a collective prop- 
erty of the potential energy surface. 

We showed that the organization of the shear trans- 
formations into cascades during the individual plastic 
events caused a scaling of the average event size with the 
length of the simulation cell regardless of the underlying 
interaction potential. These strongly correlated events, 
which involve lengths as large as the simulation cell, are 
in qualitative agreement with various numerical simula- 
tions, both atomistic [3 and mesoscopic |2^ l25l| . 
They induce scalings with the length of the system for 
various measures of event size (energy drop, stress drop, 
and elastic segment length) which would have been in- 



correctly predicted from an underlying picture of uncor- 
related localized events. Furthermore, various interac- 
tion potentials (Lennard- Jones, harmonic and hertzian 
springs) were shown to exhibit nearly identical behav- 
ior upon appropriately adjusting the energy scale for the 
given potential. 

Although, we were able to measure the properties of 
the nucleating modes quite precisely (at least in mod- 
erately sized systems), once cascades were initiated, the 
picture became quite complicated with simultaneous and 
overlapping shear transformations. In general, our work 
highlights the fundamental difficulties involved in decom- 
posing any cascade into constituent elementary events. 
Thus, important questions which are relevant to the con- 
struction of plasticity theories, such as the spatial struc- 
ture of a complete, isolated shear transformation p^ls^ . 
remain open. Even so, we were able to demonstrate that 
the predominant activity during the cascades was qual- 
itatively characteristic of local shear transformations, 
with incremental displacements and associated mechan- 
ical fields having predominantly quadrupolar character 
and the cumulative displacements and mechanical fields 
being reminiscent of narrow lines of slip. 

One of the most important directions for future study 
will be extending to finite temperatures and strain 
rates. The scaling analysis performed by Yamamoto 
and Onuki 0, [H^ showed that in sheared supercooled 
liquids, the dynamical correlation length becomes tem- 
perature independent at small enough temperatures — 
the "strong-shear" regime — whence it scales like the 
strain rate to some negative power. The recent athermal 
mesoscale model of plasticity proposed by Picard and co- 
workers j22| also exhibits a diverging length, but with a 
stronger divergence than that reported in the 2D sim- 
ulations by of Yamamoto and Onuki. Recent work by 
Langer, Liu and co-workers and Berthier and Barrat has 
suggested that imposed strain rate induces an effective 
temperature in athermal systems |88l l89l | , much the same 
as the effective temperature in unsheared supercooled liq- 
uids [90I I91I l9^ defined in terms of effective fluctuation 
dissipation relations |93j| , and consistent with the picture 
of Yamamoto and Onuki of decreasing correlations upon 
increasing either strain rate or temperature. 

This raises important questions regarding the cor- 
related motions in both the atomistic and mesoscopic 
simulations. What role do shear transformations play at 
finite temperature and strain rate both in the strong- 
shear and thermal (aging) regimes, and how robust with 
respect to finite drive is the cascade mechanism? Is 
the underlying mechanism for the heterogeneity in the 
strong-shear regime different than the thermal regime? 
Does the strong-shear scaling relation, ^ ~ 7^^^"', of 
Yamamoto and Onuki continue to hold below the glass 
transition, or is there a crossover to the ^ ^ 
behavior observed in the mesoscale model of Picard et. 
al.l More generally, is there any difference between a 
supercooled liquid in the strong-shear regime and an 
amorphous solid? These are some of the fundamental 
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questions which must be addressed in order to make 
progress' toward a coherent theory of sheared amorphous 
material. 
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